#if !defined( SPARSE_MATRIX_IO_H)
#define  SPARSE_MATRIX_IO_H
#include <iostream>
#include <vector>
#include <algorithm>
#include <cassert>
#include <string>
#include <stdexcept>
// boost headers
#include <boost/lexical_cast.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>

#include "matrix_formats.h"
#include "index_string.h"
#include "matrix.h"
#include "ultoa.h"
#include "progress_indicator.h"
#include "matrix.h"
#include "stlhelper.h"
#include "perlutil.h"
#include "progress_indicator.h"
#include "binary_io.h"
#include "tab_delimited_str.h"

#include "matrix_io.h"



class t_progress_indicator;
namespace bcpp_matrix
{


//________________________________________________________________________________________


// for reading in matrices from stream


//________________________________________________________________________________________



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	Read and write sparse matrix

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

// writing matrices

//________________________________________________________________________________________
template <typename T> 
void write_matrix(	std::ostream&				os, 
					const t_sparse_matrix<T>&	io, 
					eMATRIX_IOFORMAT			format,
					t_progress_indicator& 		dots,
					const T&					binary_null = t_matrix_traits<T>::binary_null());
// no dots
template <typename T> 
void write_matrix(	std::ostream&				os, 
					const t_sparse_matrix<T>&	io, 
					eMATRIX_IOFORMAT			format,
					const T&					binary_null = t_matrix_traits<T>::binary_null())
{
	std::ostringstream ignore;
	t_progress_indicator dots(ignore,  10000000);
	return write_matrix(os, io, format, dots, binary_null);
}

//________________________________________________________________________________________

// reading matrices

//________________________________________________________________________________________
// keeping a count of lines
template <typename T>                               		
bool read_matrix(	std::istream&				is, 
					t_sparse_matrix<T>&			io, 
					eMATRIX_IOFORMAT			format,
					unsigned&					cnt_lines,
					t_progress_indicator& 		dots,
					const T&					binary_null = t_matrix_traits<T>::binary_null());
// no cnt_lines
template <typename T>                               		
bool read_matrix(	std::istream&				is, 
					t_sparse_matrix<T>&			io, 
					eMATRIX_IOFORMAT			format,
					t_progress_indicator& 		dots,
					const T&					binary_null = t_matrix_traits<T>::binary_null())
{
	unsigned cnt_lines = 0;
	return read_matrix(is, io,  format, cnt_lines, dots, binary_null);
}

// no cnt_lines, no dots
template <typename T>                               		
bool read_matrix(	std::istream&				is, 
					t_sparse_matrix<T>&			io, 
					eMATRIX_IOFORMAT			format,
					const T&					binary_null = t_matrix_traits<T>::binary_null())
{
	std::ostringstream ignore;
	t_progress_indicator dots(ignore,  10000000);
	return read_matrix(is, io, format, dots, binary_null);
}
















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	Read and write sparse matrix


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

// input stream
template <typename T> 
struct t_sparse_matrix_istrm_wrapper
{
	t_sparse_matrix<T>&				sparse_matrix;
	eMATRIX_IOFORMAT				format;
	unsigned&						cnt_lines;
	t_progress_indicator& 			dots;
	T								binary_null;

	t_sparse_matrix_istrm_wrapper(	t_sparse_matrix<T>&		sparse_matrix_,
									eMATRIX_IOFORMAT		format_,
									unsigned&				cnt_lines_,
									t_progress_indicator& 	dots_,
									const T&				binary_null_ = 
															t_matrix_traits<T>::binary_null())
		:sparse_matrix(sparse_matrix_), format(format_), 
			cnt_lines(cnt_lines_), dots(dots_),
			binary_null(binary_null_){}
};



// output stream
template <typename T> 
struct t_sparse_matrix_ostrm_wrapper
{
	t_sparse_matrix<T>&		sparse_matrix; 
	eMATRIX_IOFORMAT		format;
	t_progress_indicator& 	dots;
	T						binary_null;

	t_sparse_matrix_ostrm_wrapper(t_sparse_matrix<T>&		sparse_matrix_,
								  eMATRIX_IOFORMAT			format_,
								  t_progress_indicator& 	dots_,
									const T&				binary_null_ = 
															t_matrix_traits<T>::binary_null())
		:sparse_matrix(sparse_matrix_), format(format_), dots(dots_),
			binary_null(binary_null_){}
};


// typedef shortcuts
typedef t_sparse_matrix_istrm_wrapper<float> t_sparse_matrix_istrm_wrapper_f;
typedef t_sparse_matrix_ostrm_wrapper<float> t_sparse_matrix_ostrm_wrapper_f;






template <typename T> 
inline std::ostream& operator <<(std::ostream& strm, const t_sparse_matrix_ostrm_wrapper<T>& wrapper)
{
	write_matrix(strm, wrapper.sparse_matrix, wrapper.format, wrapper.dots, wrapper.binary_null);
	return strm;
}
template <typename T> 
inline std::istream& operator >> (std::istream& istrm, t_sparse_matrix_istrm_wrapper<T> wrapper)
{
	if (!read_matrix(istrm,  wrapper.sparse_matrix,  wrapper.format,
					 wrapper.cnt_lines, wrapper.dots, wrapper.binary_null))
		istrm.setstate(std::ios::failbit);
	else
		// in case set at eof
		istrm.clear(istrm.rdstate() & ~(std::istream::failbit));
	return istrm;
}
















namespace matrix_io_detail
{

//________________________________________________________________________________________

//	pairs_to_matrix

//________________________________________________________________________________________
template<typename T>
void pairs_to_matrix(   std::deque<t_id_pair_value<T> >&					id_pair_values,
						typename t_sparse_matrix<T>::sparse_matrix_type&	matrix,
						t_progress_indicator&								dots,
						unsigned											matrix_size,
						bool												diagonals_used)
{
	// duplicate with rows and columns reversed
	unsigned size = id_pair_values.size();
	for (unsigned i = 0; i < size; ++i)
		id_pair_values.push_back(id_pair_values[i].reversed());
	if (!diagonals_used)
	{
		T zero_value = t_matrix_traits<T>::zero();
		for (unsigned i = 0; i < matrix_size; ++i)
			id_pair_values.push_back(t_id_pair_value<T>(zero_value, i, i));
	}
	sort(id_pair_values.begin(), id_pair_values.end());
	
	for (typename std::deque<t_id_pair_value<T> >::const_iterator i = id_pair_values.begin(),
		 end = id_pair_values.end(); i != end; ++i)
	{
		unsigned id1		= i->id_1;
		unsigned id2		= i->id_2;
		matrix(id1, id2)	= i->value;
		++dots;
	}

}

} //matrix_io_detail

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix in pairs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template<typename T>
void read_matrix_from_pairs(std::vector<std::string>::const_iterator	begin,
							std::vector<std::string>::const_iterator	end,
							t_sparse_matrix<T>&							data,
							t_progress_indicator&						dots)
{
	//
	// read pairs into array, with names in hash
	//
	std::deque<matrix_io_detail::t_id_pair_value<T> >	id_pair_values;
	t_index_string										index_names;

	// first string might be matrix name
	if (begin != end && (*begin)[0] == '>')
	{
		data.matrix_name = begin->substr(1);
		begin++;
	}
	unsigned diagonals_used = 0;
	while (begin != end)
	{
		++dots;
		diagonals_used += process_line_with_pair(*begin++, id_pair_values, index_names);
	}

	//
	//  make pairs into matrix
	//
	data.allocate(index_names.size());

	std::copy(index_names.get_deque().begin(), index_names.get_deque().end(), 
			  data.names.begin());
	matrix_io_detail::pairs_to_matrix(id_pair_values, data.matrix, dots, 
									  index_names.size(), diagonals_used);
}


template <typename T>
bool read_matrix_from_pairs(	std::istream&				infile,
								t_sparse_matrix<T>&			sparse_matrix,
								unsigned&					cnt_lines,
								t_progress_indicator&		dots)
{
	std::deque<matrix_io_detail::t_id_pair_value<T> >	id_pair_values;
	t_index_string				index_names;
	std::string line;
	//dots.set_every(20000);
	bool matrix_read = false;
	unsigned diagonals_used = 0;
	while (infile)
	{
		++cnt_lines;
		++dots;

		//
		//	get matrix name if first line
		//
		if (infile.peek() == '>')
		{
			if (!index_names.size())
			{
				infile.get();
				getline(infile, sparse_matrix.matrix_name);
				matrix_read = true;
			}

			//
			//	else matrix name belongs to next matrix
			//
			else
				break;
		}

		if(!getline(infile, line))
			break;

		if (!line.length())
		{
			// ignore consecutive blank lines
			if (!index_names.size())
				continue;

			// blank lines separate clusters
			else
				break;
		}


		diagonals_used += process_line_with_pair(line, id_pair_values, index_names);
	}

	//std::cerr << id_pair_values.size() << " pairs\n";
	//std::cerr << index_names.size() << " columns\n";
	
	if (index_names.size())
		matrix_read = true;

	if (!matrix_read)
		return false;

	//
    //  make pairs into matrix
	//
	sparse_matrix.allocate(index_names.size());

	std::copy(index_names.get_deque().begin(), index_names.get_deque().end(), 
			  sparse_matrix.names.begin());
	matrix_io_detail::pairs_to_matrix(id_pair_values, sparse_matrix.matrix, dots, 
									  index_names.size(), diagonals_used);

	return true;
}  /* read_matrix_from_stream */



template<typename T>
void write_matrix_to_pairs(std::ostream& os, const t_sparse_matrix<T>& sparse_matrix, t_progress_indicator& dots)
{
	assert(sparse_matrix.size() == sparse_matrix.names.size());

	// dots
	//dots.use_timer();
	//dots.set_limit(sparse_matrix.get_count_pairs());
	//dots.set_every(sparse_matrix.get_count_pairs() / 50);

	if (sparse_matrix.matrix_name.length())
		os << ">" << sparse_matrix.matrix_name << "\n";


	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator1 i1_t;
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator2 i2_t;

	// print out sparse rows
	for (i1_t i1 = sparse_matrix.matrix.begin1(), 
		 end1 = sparse_matrix.matrix.end1(); i1 != end1; ++i1)
		for (i2_t i2 = i1.begin(); i2 != i1.end(); ++i2)
		{
			++dots;
			os	<< sparse_matrix.names[i2.index2()] << "\t"  
				<< sparse_matrix.names[i2.index1()] << "\t" 
				<< *i2 << "\n";
		}
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix in binary

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template <typename T>
void write_matrix_to_binary(std::ostream& os, const t_sparse_matrix<T>& sparse_matrix, 
									t_progress_indicator& dots,
									T binary_null)
{
	assert(sparse_matrix.size() == sparse_matrix.names.size());


	binary_io::binary_write<unsigned>(os, 12345678);

	binary_io::binary_write(os, sparse_matrix.matrix_name);
	binary_io::binary_write(os, sparse_matrix.size());

	// use pairwise format if sparse matrix is very sparse
	bool pair_wise_format = sparse_matrix.get_count_missing() * 1.0 / 
									sparse_matrix.get_count_pairs() > 0.5;
	binary_io::binary_write<unsigned>(os, pair_wise_format ? 1 :0);

	for (unsigned i = 0, sz = sparse_matrix.size(); i < sz; ++i)
		binary_io::binary_write(os, sparse_matrix.names[i]);

	

	// 	iterator types for transversal
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator1 t_i1;
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator2 t_i2;


	// debug check
	unsigned cnt_written = 0;

	// data is very sparse: use pairwise format
	if (pair_wise_format)
	{

		// write preamble
		binary_io::binary_write(os, sparse_matrix.get_count_filled());

		
		t_i1 i1 = sparse_matrix.matrix.begin1();
		t_i1 end1 = sparse_matrix.matrix.end1();
		for (;i1 != end1; ++i1)
			for (t_i2 i2 = i1.begin(); i2 != i1.end(); ++i2)
			{
				++dots;
				++cnt_written;
				binary_io::binary_write<unsigned short>(os, i2.index1());
				binary_io::binary_write<unsigned short>(os, i2.index2());
				binary_io::binary_write<T>(os, *i2);
			}
	}


	// data is not very sparse: use dense format
	else
	{
		unsigned matrix_size = sparse_matrix.size();

		// last row index written
		unsigned last_ii = 0;

		t_i1 i1		= sparse_matrix.matrix.begin1();
		t_i1 end1	= sparse_matrix.matrix.end1();
		for (;i1 != sparse_matrix.matrix.end1(); ++i1)
		{	
			//	write empty rows
			for (unsigned jj = 0; jj < matrix_size * (last_ii - i1.index1()); ++jj)
				binary_io::binary_write<T>(os, binary_null);
			last_ii = i1.index1() + 1;

			// last column index written
			unsigned last_jj = 0;
			
			for (t_i2 i2 = i1.begin(); i2 != i1.end(); ++i2)
			{
				//	write empty columns
				for (; last_jj != i2.index2(); ++last_jj)
					binary_io::binary_write<T>(os, binary_null);

				// write non-empty columns
				++dots;
				++cnt_written;
				++last_jj;
				binary_io::binary_write<T>(os, *i2);
			}

			//	write empty columns
			for (; last_jj != matrix_size; ++last_jj)
				binary_io::binary_write<T>(os, binary_null);
		}
	}
	assert(cnt_written == sparse_matrix.get_count_filled());

}



template <typename T>
bool read_matrix_from_binary(std::istream& 				is, 		
							t_sparse_matrix<T>&			sparse_matrix,
							unsigned&					cnt_lines,
							t_progress_indicator&		dots,
							T							binary_null)
{

	// ignore blank lines at start
	while (is.peek() == '\n')
	{
		++cnt_lines;
		is.ignore(1);
	}
	if (is.eof())
		return false;

	unsigned check = binary_io::binary_read<unsigned>(is);
	if (check != 12345678)
		throw std::runtime_error("Matrix in binary format is corrupted and cannot be read.\n");

	sparse_matrix.matrix_name = binary_io::binary_read<std::string>(is);

	unsigned sz = binary_io::binary_read<unsigned>(is);
	sparse_matrix.allocate(sz);
	bool pair_wise_format = binary_io::binary_read<unsigned>(is) != 0;

	for (unsigned i = 0; i < sz; ++i)
		sparse_matrix.names[i] = binary_io::binary_read<std::string>(is);


	if (pair_wise_format)
	{
		unsigned cnt_read = binary_io::binary_read<unsigned>(is);
		while (cnt_read)
		{
			unsigned i = binary_io::binary_read<unsigned short>(is);
			unsigned j = binary_io::binary_read<unsigned short>(is);
			sparse_matrix.matrix(i, j) = binary_io::binary_read<T>(is);
			--cnt_read;
			++dots;
		}
	}

	// dense format
	else
	{
		for (unsigned i = 0; i < sz; ++i)
		{
			++dots;
			for (unsigned j = 0; j < sz; ++j)
			{
				T val = binary_io::binary_read<T>(is);
				if (val != binary_null)
					sparse_matrix.matrix(i, j) = val;
			}
		}
	}
	return true;
}










//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix in square

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
			
//________________________________________________________________________________________
//
//  write to square
// 
//________________________________________________________________________________________
template<typename T>
void write_matrix_to_square(std::ostream&					os, 
							const t_sparse_matrix<T>&		sparse_matrix, 
							bcpp_matrix::eMATRIX_IOFORMAT	format, 
							t_progress_indicator&			dots,
							T								binary_null)
{
	assert(sparse_matrix.size() == sparse_matrix.names.size());

	// no preamble for phylip
	if (format != bcpp_matrix::eMIO_PHYLIP && format != bcpp_matrix::eMIO_PHYLIP_NULL)
	{
		if (sparse_matrix.matrix_name.length())
			os << ">" << sparse_matrix.matrix_name << "\n";
	}

	os << sparse_matrix.size() << '\n';


	std::string str_null = '\t' + get_null_str(binary_null);
	std::string non_null;

	// add null weighting
	if (format == bcpp_matrix::eMIO_PHYLIP_NULL)
	{
		str_null +="\t0";
		non_null +="\t1";
	}

	unsigned matrix_size = sparse_matrix.size();

	unsigned cnt_written = 0;

	// 	iterator types for transversal
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator1 t_i1;
	typedef typename t_sparse_matrix<T>::sparse_matrix_type::const_iterator2 t_i2;
	t_i1 i1		= sparse_matrix.matrix.begin1();
	t_i1 end1	= sparse_matrix.matrix.end1();

	// last row index written
	unsigned last_ii = 0;

	for (;i1 != end1; ++i1)
	{	
        //
        //  write empty rows
        //
        unsigned non_empty_row = i1.index1();
		for (; last_ii != non_empty_row; ++last_ii)
		{
			// write name followed by missing flags
			os << phylip_pad_name(sparse_matrix.names[last_ii], format);
			// lower doesn't include diagonals
			unsigned j_sz = format == bcpp_matrix::eMIO_LOWER  ? last_ii : matrix_size;
			for (unsigned j = 0; j < j_sz; ++j)
				os  << str_null;
			os << '\n';
		}
		
		// write name
		os << phylip_pad_name(sparse_matrix.names[non_empty_row], format);


		// last column index written
		unsigned last_jj = 0;

		unsigned j_sz = format == bcpp_matrix::eMIO_LOWER  ? last_ii : matrix_size;
		for (t_i2 i2 = i1.begin(); i2 != i1.end(); ++i2)
		{
			if (i2.index2() >= j_sz)
				break;
			
			//	write empty columns
			for (; last_jj != i2.index2(); ++last_jj)
				os  << str_null;

			// write non-empty column
			++dots;
			os << '\t' << *i2 << non_null;
			++cnt_written;
			++last_jj;

		}

		//	write empty columns
		for (; last_jj != j_sz; ++last_jj)
			os  << str_null;


		// update last ii for writing non_empty_row
		++last_ii;
		os << '\n';
	}
	//assert(cnt_written == sparse_matrix.get_count_filled());


}

//________________________________________________________________________________________
//
//  read from square
// 
//________________________________________________________________________________________
template <typename T>
bool read_matrix_from_square(	std::istream&					is, 
								t_sparse_matrix<T>&				sparse_matrix,
								bcpp_matrix::eMATRIX_IOFORMAT	format, 
								unsigned&						cnt_lines,
								t_progress_indicator&			dots,
								const T&						binary_null)
{

	// ignore blank lines at start
	while (is.peek() == '\n')
	{
		++cnt_lines;
		is.ignore(1);
	}
	if (is.eof())
		return false;

	bool matrix_read = false;
	if (is.peek() == '>')
	{
		is.ignore(1);
		++cnt_lines;
		getline(is, sparse_matrix.matrix_name);
		matrix_read = true;
	}

	if (is.eof())
		return matrix_read;

	unsigned sz;
	is >> sz;
	if (!is)
		throw std::runtime_error("Invalid square matrix format.\nThe matrix size on the first line is missing.");
	++cnt_lines;
	sparse_matrix.allocate(sz);
	std::string str_null = get_null_str(binary_null);
	for (unsigned i = 0; i < sz; ++i)
	{
		++dots;
		char endc;
		is >> tab_delimited_str('\n', sparse_matrix.names[i], endc);
		++cnt_lines;
		unsigned j_sz = format == bcpp_matrix::eMIO_LOWER  ? i : sz;
		for (unsigned j = 0; j < j_sz; ++j)
		{
			// read first identifier
			std::string value;
			char endc;
			is >> tab_delimited_str('\t', value, endc);
			if (value == "")
				throw std::runtime_error("Invalid square matrix_format.\nEmpty value on line " +
   											ultoa(cnt_lines));
			// phylip replicates
			if (format == bcpp_matrix::eMIO_PHYLIP_NULL)
			{
				std::string replicate_value;
				is >> tab_delimited_str('\t', replicate_value, endc);
				if (replicate_value == "")
					throw std::runtime_error("Invalid square matrix_format.\nEmpty replicate number on line " +
												ultoa(cnt_lines));
				int cnt_replicates = t_matrix_traits<int>::from_string(replicate_value);
				if (!cnt_replicates)
					continue;
			}

			// only save if not null
			if (value != str_null)
			{
				sparse_matrix.matrix(i, j) = t_matrix_traits<T>::from_string(value);

				// lower format is the only one which is symmetrical
				if (format == bcpp_matrix::eMIO_LOWER)
					sparse_matrix.matrix(j, i) = sparse_matrix.matrix(i, j);
			}

		}
		if (!is)
			throw std::runtime_error("Invalid square matrix format.");
	}
	T zero_value = static_cast<T>(0U);
	if (format == bcpp_matrix::eMIO_LOWER)
		for (unsigned i = 0; i < sparse_matrix.size(); ++i)
			sparse_matrix.matrix(i, i) = zero_value;
	return true;
}












//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	all formats

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template<typename T>
void write_matrix(	std::ostream&				os, 
					const t_sparse_matrix<T>&	sparse_matrix,
					eMATRIX_IOFORMAT			format,
					t_progress_indicator& 		dots,
					const T&					binary_null)
{
	switch (format)
	{
		case eMIO_BINARY:
			write_matrix_to_binary(os, sparse_matrix, dots, binary_null);
			break;
		case eMIO_PAIRS :
			write_matrix_to_pairs(os, sparse_matrix, dots);
			break;
		case eMIO_LOWER :
		case eMIO_SQUARE:
		case eMIO_PHYLIP:
		case eMIO_PHYLIP_NULL:
			write_matrix_to_square(os, sparse_matrix, format, dots, binary_null);
			break;
		default:
			assert(1);
	}
}

template <typename T>
bool read_matrix(	std::istream&				is,
					t_sparse_matrix<T>&			sparse_matrix,
					eMATRIX_IOFORMAT			format,
					unsigned&					cnt_lines,
					t_progress_indicator& 		dots,
					const T&					binary_null)
{
	try
	{
		switch (format)
		{
			case eMIO_BINARY:
				return read_matrix_from_binary(is, sparse_matrix, cnt_lines, dots, binary_null);
			case eMIO_PAIRS :
				return read_matrix_from_pairs(is, sparse_matrix, cnt_lines, dots);
			case eMIO_LOWER :
				return read_matrix_from_square(is, sparse_matrix, format, cnt_lines, dots, binary_null);
			case eMIO_SQUARE:
			case eMIO_PHYLIP:
			case eMIO_PHYLIP_NULL:
				return read_matrix_from_square(is, sparse_matrix, format, cnt_lines, dots, binary_null);
			default:
				assert(1);
		}
	}
	catch (std::runtime_error& error)
	{
		throw std::runtime_error(std::string(error.what()) + "\nThe error occurred at line # " + ultoa(cnt_lines + 1));
	}

	throw std::logic_error("Write matrix: unknown format\n");
	return true;
}

};



#endif	//SPARSE_MATRIX_IO_H

